
function [liquidCompPost,augite, temperature,phaseProportions_out,PhaseChangeInfo,liquidCompPost_Trace,solidphases,DistanceMatrix,OPAMout,OPALMout,olivineKD,PO_Tormey_liquids_out,Fout] = crystallizeFUNCTION_while4_fixed(liqCompIn, P,Pressures_OP,percentIncrement,targetCPXFoContent, targetPlagFoContent, targetFoContent,errorOPC, errorOP,magnetiteInfo, whichRegression,distanceValue,Elements,Weights, Components4OPAM, TemperatureIndicies, CompIndicies, PressureIndicies, FittedVariablesIndicies, ModelVariablesIndicies,ElementIndicies4Target_Tormey,ModelVariablesIndiciesOPLAM,ElementIndicies4Target_Yang,ElementIndicies4Target_OUT,Trace,Dmatrix)

olivineKD = 0.30; % is updated later
MASTER_REGRESSIONS
AcceptableTotalRange_HERE=[98.5 101.5]; 

% FailedVariable = 'GOOD';
if nargin < 25
    Dmatrix=[1 1 1 1 1 1 1];
    liquidCompPost_Trace=[];
end

% outputs:
% %['1-SiO2'	'2-TiO2' '3-Al2O3'	'4-Cr2O3'  '5-FeO'	'6-MnO'	'7-MgO'	'8-CaO'  '9-Na2O'	'10-K2O' '11-P2O5' '12-SUM 13-MG#MELT 14-MG#OLIV 15-NAK#MELT 16-FC_EXTENT]

% the possible phases to crystallize or melt are:
% [olivine opx pigeonite augite plagioclase garnet amphibole ilmenite magnetite]


%               1       2       3         4      5        6     7      8          9         10         11
targetElements4minerals ={ 'SiO2'	'TiO2' 'Al2O3'	'Cr2O3'  'FeO'	'MnO'	'MgO'	'CaO'    'Na2O'	'K2O'	    'P2O5'};
[A,ElementIndicies4Minerals] = ismember(targetElements4minerals, Elements);



targetElements_Tormey ={ 'SiO2'	'TiO2' 'Al2O3'	'Cr2O3'  'FeO'	'MnO'	'MgO'	'CaO'    'Na2O'	'K2O' 'P2O5' 'Fe2O3'};
[A,ElementIndicies4Minerals_Tormey] = ismember(targetElements_Tormey, targetElements4minerals);

[A,ElementIndicies4Target_Sum] = ismember(Elements, targetElements4minerals);


[A,MgNumCol] = ismember('Mg#', Elements);
[A,NaKNumCol] = ismember('NaK#', Elements);



%% STEP 1: change inputs/outputs for your needs.
%%%  For example, if you want to save the olivine composition too, just change the fucntion to
%%% [liquidCompPost,F,olivine] = crystallizeFUNCTION(liqComp, phaseProportions, percent,percentIncrement)



%% STEP 2: check KDs
% For now we are using these set KDs. Change them to your liking.
% In the future we could have a parameterized value that calculates within
% the routine based on the liquid comp (or T and P).
%%%     oliv   opx  pig/augite  plag  amphibole garnet[FeMg  SiMg]
%%%      1      2           3         4          5                   6       7
KDs = [olivineKD  0.27      0.23      1       0.35             0.48   0.8];

%% STEP 3: define any composition-invariant phases
%            1       2       3         4      5        6     7      8          9         10         11
%          SiO2	TiO2 Al2O3	Cr2O3  FeO	MnO	MgO	CaO	    Na2O	K2O	    P2O5
%%magnetite and ilmenite comps are already defined.  You can set to known comp for any
%%mineal that is not included above
%magnetiteComp = [0   7.8  2.7 0   87.9  0   1.6  0   0  0   0]; %%This is the default in the fortran code, change this for your comp!
MagnetiteComp = [0.34	18.5	2.49	0	67.9	0.49	3.9	0	0	0	0];
ilmeniteComp = [0 45 0.5 0 50 0.5 4 0 0 0 0]; %%This is some real comp from Ben's expts, change!
blankComp = [0 0 0 0 0 0 0 0 0 0 0];
augite=blankComp;
temperature=[];
%%% pure ilmenite FeTiO3 = [0 52.65 0 0 47.35 0 0 0 0 0 0];

%% Now calculate crystallization:

%%% the liquid composition matrix initially saves the original liq comp in
%%% the first row, edit line 60 (**) to change this
%liquidCompPost(1,:) = liqComp; %% (**): change here if you don't want to save the initial liq in the derivative liquid comp matrix
n =1; % i and n will be different if the percentIncrmenet is not = 1
Increment = percentIncrement/100;
F(1) = 1;


% liqCompIn
% liqCompIn(ElementIndicies4Minerals) = liqCompIn(ElementIndicies4Minerals)./nansum(liqCompIn(ElementIndicies4Minerals))*100; 
%nansum(liqCompIn(ElementIndicies4Minerals))


liquidCompPost(1,:) = liqCompIn;

liquidCompPost_Trace(1,:) = Trace;
phaseProportions_out(1,:)  = NaN.*[0       0         0         0        0             0        0            0              0];

InitialLiquid_Tormey = TormeyProjection(liquidCompPost(1,:),Elements,ElementIndicies4Target_Tormey);%[ 1-Quartz 2-Plag 3-Olivine 4-Cpx]
PlagOliv_InitialLiquid = PlagOlivProps_QPOC(InitialLiquid_Tormey);
PO_Tormey_liquids_out(1,:) = PlagOliv_InitialLiquid;
OPAM_initialLiquid_RMSD=[];

TormeyOPAM_MEcomp= OPAM_Yang1996(liquidCompPost(1,:),P,Elements,ElementIndicies4Target_Yang,ElementIndicies4Target_OUT);
TormeyOPAM= TormeyProjection(TormeyOPAM_MEcomp, Elements,ElementIndicies4Target_Tormey);
PlagOliv_OPAM = PlagOlivProps_QPOC(TormeyOPAM);

% Distance from OPAM
OPAM_initialLiquid_RMSD= RMSD(InitialLiquid_Tormey, TormeyOPAM, Components4OPAM, 'NormalizeComponents');
% Distance from OP-OPAM
OP_OPAM_initialLiquid_RMSD= RMSD(PlagOliv_InitialLiquid, PlagOliv_OPAM, [1 2], 'NormalizeComponents');


% This is calculating the OPALM point at the pressure of FC in the plag
% field in order to calculate the OP boundary by projecting the OPALM
% boundary
[OPALM] = calculate_OPALM_MASTER(liquidCompPost(1,:),P,TemperatureIndicies,CompIndicies,ModelVariablesIndiciesOPLAM,...
    eval(sprintf('garnet_coefficients_%s',whichRegression)), eval(sprintf('spinel_coefficients_%s',whichRegression)), eval(sprintf('plagioclase_coefficients_%s',whichRegression)));


%         [OPALM_allP] = calculate_OPALM_MASTER(liquidCompPost(1,:),Pressures_OP,Elements,...
%                     eval(sprintf('VariablesT_%s',whichRegression)), eval(sprintf('VariablesX_%s',whichRegression)), eval(sprintf('VariablesP_%s',whichRegression)),...
%                     eval(sprintf('garnet_coefficients_%s',whichRegression)), eval(sprintf('spinel_coefficients_%s',whichRegression)), eval(sprintf('plagioclase_coefficients_%s',whichRegression)));

test = XPTmeter_MASTER(liquidCompPost(1,:),TemperatureIndicies, CompIndicies, PressureIndicies, FittedVariablesIndicies, ModelVariablesIndicies,...
    eval(sprintf('garnet_coefficients_%s',whichRegression)), eval(sprintf('spinel_coefficients_%s',whichRegression)), eval(sprintf('plagioclase_coefficients_%s',whichRegression)));


OPALM_Melting_Distance = test.fits;
OPALM_Melting_Distance_PT = test.PTTmp(1:2);
OPALM_Melting_Distance_P = test.PTTmp([1 4 7]);
OPALM_Melting_Distance_T = test.PTTmp([2 5 8]);
OPALM = OPALM(3:end);

PlagOliv_OPALM = PlagOlivProps_QPOC(OPALM);
%PO_OPALM_allP = PlagOlivProps_QPOC(OPALM_allP(:,3:end));


% Distance from OPALM
OPALM_initialLiquid_RMSD= RMSD(InitialLiquid_Tormey, OPALM, Components4OPAM, 'NormalizeComponents');
% Distance from OP-OPALM
OP_OPALM_initialLiquid_RMSD= RMSD(PlagOliv_InitialLiquid, PlagOliv_OPALM, [1 2], 'NormalizeComponents');

Garnet2SpinelTransition = (1.2-0.52.*liquidCompPost(NaKNumCol)+2.11.*liquidCompPost(MgNumCol))*10;
Spinel2PlagTransition = (0.823+1.139.*liquidCompPost(NaKNumCol)+0.165.*liquidCompPost(MgNumCol))*10;


DistanceMatrix(1,:) = [OPAM_initialLiquid_RMSD OP_OPAM_initialLiquid_RMSD OPALM_initialLiquid_RMSD OP_OPALM_initialLiquid_RMSD OPALM_Melting_Distance ...
    InitialLiquid_Tormey(4)-TormeyOPAM(4) Garnet2SpinelTransition Spinel2PlagTransition OPALM_Melting_Distance_P OPALM_Melting_Distance_T];


OPAMout(1,:) = TormeyOPAM;
OPALMout(1,:) = OPALM;

% PlagOliv_OPAM
% PlagOliv_OPALM
% PlagOliv_InitialLiquid
% PO_OPALM_allP


while  liquidCompPost(n,MgNumCol) < targetFoContent
    
    
    [olivine(n,:), olivineKD(n,:)]= fractionateOlivine_V2(liquidCompPost(n,ElementIndicies4Minerals),KDs(1));
    plagioclase(n,:) = fractionatePlagioclase_V2(liquidCompPost(n,ElementIndicies4Minerals),KDs(4));
    
    %GROVE MORBFRAC
    % augite(n,:) = fractionateAugite_V2(liquidCompPost(n,1:numElements),targetElements, P,KDs(3));
    
    % XX YANG CPX
    [augite(n,:),temperature(n,:)] = fractionateAugite_V3_Yang1996(liquidCompPost(n,ElementIndicies4Minerals), P,KDs(3));
    
    
    %SMB: save time by not calcualting these compositions. Remember to uncomment if you want to use them!
    %     opx(n,:) = fractionateOpx_V2(liquidCompPost(n,1:numElements),KDs(2));
    %     pigeonite(n,:) = fractionatePIG_V2(liquidCompPost(n,1:numElements),KDs(3));
    %     amphibole(n,:) = fractionateAmphibole_V2(liquidCompPost(n,1:numElements),KDs(5));
    %     garnet(n,:) = fractionateGarnet_V2(liquidCompPost(n,1:numElements),KDs(6),KDs(7));
    %     %%% uses input compositions, can change later to make it liquid
    %     %%% comp, T,P, fO2 dependent if useful
    %     ilmenite(n,:) = ilmeniteComp;
    %     magnetite(n,:) = magnetiteComp;
    opx(n,:) = blankComp;
    pigeonite(n,:) = blankComp;
    amphibole(n,:) = blankComp;
    garnet(n,:) =blankComp;
    ilmenite(n,:) =blankComp;
    magnetite(n,:) = MagnetiteComp;
    
    %%% now combine the mineral data into a matrix
    %%% each row is a mineral composition.  If you add minerals,
    %%% make sure to change it here.
    %%% [olivine opx  pigeonite augite plagioclase garnet amphibole ilmenite magnetite]
    
    %clearvars temp
    
    temp = [olivine(n,:); opx(n,:); pigeonite(n,:); augite(n,:); plagioclase(n,:); garnet(n,:); amphibole(n,:); ilmenite(n,:); magnetite(n,:)];
    
    
    TormeyMineralComps=[olivine(n,:); plagioclase(n,:); augite(n,:)];
    
    
    TormeyMineralComps_actaullyTormey= TormeyProjection(TormeyMineralComps,targetElements4minerals,ElementIndicies4Minerals_Tormey);
    TormeyLiquid_actaullyTormey= TormeyProjection(liquidCompPost(n,:),Elements,ElementIndicies4Target_Tormey);
    
    TormeyMineralComps_actaullyTormey(:,[1 5:7]) = [];
    TormeyLiquid_actaullyTormey(:,[1 5:7])  = [];
    
    TormeyMineralComps_actaullyTormey = bsxfun(@times, TormeyMineralComps_actaullyTormey,nansum(TormeyMineralComps_actaullyTormey,2));
    TormeyLiquid_actaullyTormey = bsxfun(@times, TormeyLiquid_actaullyTormey,nansum(TormeyLiquid_actaullyTormey,2));
    
    [MassBalPP] = (TormeyMineralComps_actaullyTormey'\TormeyLiquid_actaullyTormey');
    MassBalPP=MassBalPP./sum(MassBalPP);
    
    % XX MASS BALANCE PP Tormey Components, could change to oxides
    phaseProportions =  [-MassBalPP(1)       0         0         -MassBalPP(3)       -MassBalPP(2)             0        0            0              0];
    % mass conversion???
    
    distanceMatrixValue2Check = 1;
    distanceMatrixValue2Check_error = errorOPC;
    
    if liquidCompPost(n,MgNumCol) > targetCPXFoContent %Switch to OP!
        % 'Qtz'    'Plag'    'Oliv'    'Cpx'
        TormeyLiquid_actaullyTormey= TormeyProjection(liquidCompPost(n,:),Elements,ElementIndicies4Target_Tormey);
        PlagOl_Liquid = PlagOlivProps_QPOC(TormeyLiquid_actaullyTormey);
        
        PlagOl_Liquid = [PlagOl_Liquid(1).*2.7 PlagOl_Liquid(2).*3.2];
        PlagOl_Liquid = PlagOl_Liquid./sum(PlagOl_Liquid);
        
        
        %USE THE Liquid COORDIANTES
        OLPROP = PlagOl_Liquid(2);
        PLPROP=  PlagOl_Liquid(1);
        phaseProportions =  [-1*OLPROP       0         0        0      -1*PLPROP             0        0            0              0];
        augite(n,:)=0.*augite(n,:);
        distanceMatrixValue2Check = distanceValue;
        distanceMatrixValue2Check_error = errorOP;
    end
    
    
    if liquidCompPost(n,MgNumCol) > targetPlagFoContent % Switch to O-only
        %if abs(DistanceMatrix(n,distanceValue)) > OP_distance
        phaseProportions  = [-1       0         0         0        0             0        0            0              0];
        augite(n,:)=0.*augite(n,:);
        distanceMatrixValue2Check = 0;
        distanceMatrixValue2Check_error=0;
    end
    
    
    
    phaseProportions=-1.*phaseProportions./sum(phaseProportions);
    
    if liquidCompPost(n,7) < magnetiteInfo(1)
        phaseProportions(end) = -1*magnetiteInfo(2);
        phaseProportions(1:end-1) = phaseProportions(1:end-1)./(1-phaseProportions(end));
    end
    
    
    phaseProportions_out(n+1,:)=phaseProportions;
    %      TormeyLiquid = [POLIV,PPLAG,PCPX];
    %      TormeyOPAMComp = [PHJOLIV, PHJPLAG,PHJCPX];
    %      TormeyMineralComps = [POlOLIV, POlPLAG,POlCPX; ...
    %                                        PPlagOLIV, PPlagPLAG PPlagCPX;...
    %                                        PCPXOLIV, PCPXPLAG,PCPXCPX];
    %
    % 'Mass Balance: Mineral comps \Liquid comp OL PLAG CPX'
    % %calcualtes phase proportions using liquid composition [ol plag cpx]
    % FractionatingModeHere2 = (TormeyMineralComps'\TormeyLiquid') % 3x3' \ 1x3' = 3x3\3x1 = 3x1
    
    solidphases(n,:) = phaseProportions*temp; %%multiples each mineral times its proportion, saves it in a new matrix
    
    %%% the new liquid compostiion is calculated by removing an increment of
    %%% solid material.
    
    
  % phaseProportions
 %   tempLiquid = (liquidCompPost(n,ElementIndicies4Minerals) - Increment*solidphases(n,:))./(1+Increment);
     tempLiquid = (liquidCompPost(n,ElementIndicies4Minerals) - Increment*solidphases(n,:));   
% 
%     'solids'
%     solidphases(n,:)
%     nansum(solidphases(n,:))
%  
%     
%     'liquid'
%     liquidCompPost(n,ElementIndicies4Minerals) 
%     nansum(liquidCompPost(n,ElementIndicies4Minerals) )
%     
%     'new liquid'
%     tempLiquid
%     nansum(tempLiquid)
%     
%     tempLiquid = 100.*tempLiquid./nansum(tempLiquid)
%    
    
  %  return
    
    liquidCompPost(n+1,:) = sumMgNumNorm_PT(tempLiquid,targetElements4minerals,ElementIndicies4Target_Sum,AcceptableTotalRange_HERE);
    F(n+1) = (1-Increment)^n;
    
    
    
    %temp = [olivine(n,:); opx(n,:); pigeonite(n,:); augite(n,:); plagioclase(n,:); garnet(n,:); amphibole(n,:); ilmenite(n,:); magnetite(n,:)];
    
    
    % These need to be in the order of DMatrid :CPX	OPX	OLIVINE	PLAG
    % SPINEL	GARNET	Magnetite
    
    
    TracePhaseProportions=[phaseProportions(4) phaseProportions(2) phaseProportions(1) phaseProportions(5) 0 phaseProportions(6) phaseProportions(9)];
    phaseProportionsTrace_out(n,:)=TracePhaseProportions;
    D = TracePhaseProportions*Dmatrix';
    
    
    liquidCompPost_Trace(n+1,:) = (liquidCompPost_Trace(n,:) - Increment.*(liquidCompPost_Trace(n,:).*D))./(1+Increment) ;
    
    InitialLiquid_Tormey = TormeyProjection(liquidCompPost(n+1,:),Elements,ElementIndicies4Target_Tormey);   %[ 1-Quartz 2-Plag 3-Olivine 4-Cpx]
    PlagOliv_InitialLiquid = PlagOlivProps_QPOC(InitialLiquid_Tormey);
    
    
    OPAM_initialLiquid_RMSD=[];
    
    
    
    TormeyOPAM_MEcomp= OPAM_Yang1996(liquidCompPost(n+1,:),P,Elements,ElementIndicies4Target_Yang,ElementIndicies4Target_OUT);
    TormeyOPAM= TormeyProjection(TormeyOPAM_MEcomp,Elements,ElementIndicies4Target_Tormey);
    PlagOliv_TormeyOPAM = PlagOlivProps_QPOC(TormeyOPAM);
    
    
    [OPALM] = calculate_OPALM_MASTER(liquidCompPost(n+1,:),P,TemperatureIndicies,CompIndicies,ModelVariablesIndiciesOPLAM, ...
        eval(sprintf('garnet_coefficients_%s',whichRegression)), eval(sprintf('spinel_coefficients_%s',whichRegression)), eval(sprintf('plagioclase_coefficients_%s',whichRegression)));
    
    
    test = XPTmeter_MASTER(liquidCompPost(n+1,:),TemperatureIndicies, CompIndicies, PressureIndicies, FittedVariablesIndicies, ModelVariablesIndicies,...
        eval(sprintf('garnet_coefficients_%s',whichRegression)), eval(sprintf('spinel_coefficients_%s',whichRegression)), eval(sprintf('plagioclase_coefficients_%s',whichRegression)));
    OPALM_Melting_Distance = test.fits;
    
    
    OPALM_Melting_Distance_PT = test.PTTmp(1:2);
    OPALM_Melting_Distance_P = test.PTTmp([1 4 7]);
    OPALM_Melting_Distance_T = test.PTTmp([2 5 8]);
    
    OPALM = OPALM(3:end);
    OP_OPALM = PlagOlivProps_QPOC(OPALM);
    
    
    OPAM_initialLiquid_RMSD= RMSD(InitialLiquid_Tormey,TormeyOPAM,Components4OPAM,'NormalizeComponents');
    OP_OPAM_initialLiquid_RMSD= RMSD(PlagOliv_InitialLiquid, PlagOliv_TormeyOPAM, [1 2],'NormalizeComponents') ;
    
    OPALM_initialLiquid_RMSD= RMSD(InitialLiquid_Tormey,OPALM,Components4OPAM,'NormalizeComponents');
    OP_OPALM_initialLiquid_RMSD= RMSD(PlagOliv_InitialLiquid, OP_OPALM, [1 2],'NormalizeComponents') ;
    
    
    
    Garnet2SpinelTransition = (1.2-0.52.*liquidCompPost(n+1,NaKNumCol)+2.11.*liquidCompPost(n+1,MgNumCol))*10;
    Spinel2PlagTransition = (0.823+1.139.*liquidCompPost(n+1,NaKNumCol)+0.165.*liquidCompPost(n+1,MgNumCol))*10;
    
    
    
    
    OPAMout(n+1,:) = TormeyOPAM;
    OPALMout(n+1,:) = OPALM;
    DistanceMatrix(n+1,:) = ...
        [OPAM_initialLiquid_RMSD OP_OPAM_initialLiquid_RMSD OPALM_initialLiquid_RMSD OP_OPALM_initialLiquid_RMSD ...
        OPALM_Melting_Distance InitialLiquid_Tormey(4)-TormeyOPAM(4) ...
        Garnet2SpinelTransition Spinel2PlagTransition OPALM_Melting_Distance_P OPALM_Melting_Distance_T];
    
        PO_Tormey_liquids_out(n+1,:) = PlagOliv_InitialLiquid;
    
    
    %a_PO_OPALM=find(OP_OPALM_ALL_P(:,1)<PlagOliv_InitialLiquid(1) & PO_OPALM_initialLiquid_RMSD < ErOP);
    
    
    
    if distanceMatrixValue2Check>0 && DistanceMatrix(n+1,distanceMatrixValue2Check) > distanceMatrixValue2Check_error
        %         n
        %         distanceMatrixValue2Check
        %         DistanceMatrix(n+1,distanceMatrixValue2Check)
        %         distanceMatrixValue2Check_error
        %         liquidCompPost
        
        liquidCompPost = liquidCompPost.*NaN;
        augite=augite.*NaN;
        temperature=temperature.*NaN;
        
        
        liquidCompPost_Trace=liquidCompPost_Trace.*NaN;
        solidphases=solidphases.*NaN;
        %DistanceMatrix
        OPAMout=OPAMout.*NaN;
        OPALMout = OPALMout.*NaN;
        olivineKD=olivineKD.*NaN;
        
        PO_Tormey_liquids_out=PO_Tormey_liquids_out.*NaN;
        phaseProportions_out = phaseProportions_out.*NaN;
        %FailedVariable = 'FAILED';
        Fout=NaN;
        PhaseChangeInfo=NaN;
        return
    end
    
    
    
    
    
    %     if distanceMatrixValue2Check ==1  && InitialLiquid_Tormey(4)<TormeyOPAM(4)
    %
    %         'okay'
    %         n
    %         liquidCompPost
    %         DistanceMatrix
    %
    %         InitialLiquid_Tormey
    %         TormeyOPAM
    %
    %
    %         liquidCompPost = liquidCompPost.*NaN;
    %         augite=augite.*NaN;
    %         temperature=temperature.*NaN;
    %
    %
    %         liquidCompPost_Trace=liquidCompPost_Trace.*NaN;
    %         solidphases=solidphases.*NaN;
    %         %DistanceMatrix
    %         OPAMout=OPAMout.*NaN;
    %         OPALMout = OPALMout.*NaN;
    %         olivineKD=olivineKD.*NaN;
    %
    %         PO_Tormey_liquids_out=PO_Tormey_liquids_out.*NaN;
    %         phaseProportions_out = phaseProportions_out.*NaN;
    %         %FailedVariable = 'FAILED';
    %         Fout=NaN;
    %         PhaseChangeInfo=NaN;
    %         return
    %     end
    
    
    
    n=n+1;
end








% if strcmp(FailedVariable,'FAILED')==1
%
%     return
% end



Fout = [F ;(1-F)*100]';
%add on the % crystallized
%['1-SiO2'	'2-TiO2' '3-Al2O3'	'4-Cr2O3'  '5-FeO'	'6-MnO'	'7-MgO'	'8-CaO'  '9-Na2O'	'10-K2O' '11-P2O5' '12-SUM 13-MG#MELT 14-MG#OLIV 15-NAK#MELT 16-FC_EXTENT]


% [Ol]= find(phaseProportions_out(:,4)==0,1); % first step at which there is no cpx
% [OlPlag]=find(phaseProportions_out(:,4)==0 & phaseProportions_out(:,5)==0,1); % first step at which it is olivine only
% [OlPlagCpx]=find(phaseProportions_out(:,4)==0 & phaseProportions_out(:,5)==0,1); % first step at which it is olivine only



[a]= find(phaseProportions_out(:,4)==0,1); % first step at which there is no cpx
[c]=find(phaseProportions_out(:,4)==0 & phaseProportions_out(:,5)==0,1); % first step at which it is olivine only
[d]=find(phaseProportions_out(:,4)==0 & phaseProportions_out(:,5)==0,1); % first step at which it is olivine and plag only


if a>1
    PhaseChangeInfo(1) = Fout([a-1],2);
else
    PhaseChangeInfo(1) = 0;
end



if c>1
    PhaseChangeInfo(2) = Fout([c-1],2)-Fout([a],2);
    PhaseChangeInfo(3)=Fout(c,2);
else
    if a >1
        PhaseChangeInfo(2) = Fout([end-1],2)-Fout([a],2);
        PhaseChangeInfo(3) = Fout(end-1,2);
    elseif n>1
        PhaseChangeInfo(2) = Fout([end-1],2);
        PhaseChangeInfo(3) = Fout(end-1,2);
    else
        PhaseChangeInfo(2) = Fout([end],2);
        PhaseChangeInfo(3) = Fout(end,2);
    end
end

if c<=1 & n>2
    if phaseProportions_out(2,5)~=0
        PhaseChangeInfo(3) = Fout(end,2);
        PhaseChangeInfo(2) = Fout(end,2);
    end
end


if n ==1
    %DistanceMatrix(2,:)=[];
end



if n >1
    %%% converts solid phases to positive values to offset more addition than
    %%% fractionation
    solidphases = abs(solidphases);
    %%% tacks on a new column with Mg# for each compositional matrix.  Could
    %%% make this also calc NaK#s, etc...
    %liqComp= addMgNum_V2(liqComp);
    %liqComp= addMgNum_V2(olivine);
    %liquidCompPost= addMgNum_V2(liquidCompPost);
    solidphases= addMgNum_V2(solidphases);
    olivine= addMgNum_V2(olivine);
    opx= addMgNum_V2(opx);
    pigeonite= addMgNum_V2(pigeonite);
    augite= addMgNum_V2(augite);
    plagioclase= addMgNum_V2(plagioclase);
    amphibole= addMgNum_V2(amphibole);
    garnet= addMgNum_V2(garnet);
    ilmenite= addMgNum_V2(ilmenite);
    magnetite= addMgNum_V2(magnetite);
else
    
    
end
solidphases(:,MgNumCol) = (1-F(2:end))*100;
%solidphases(:,MgNumCol+1) = NaN.*solidphases(:,MgNumCol+1) ;


end

